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REAL TIME ESTIMATION AND PREDICTION 
OF SHIP MOTIONS USING KALMAN FILTER TECHNIQUES 

ABSTRACT 

A study of the real time estimation and prediction 
of ship motions, velocities and accelerations is presented. 
The ship motion estimations are of particular interest for 
operations in rough seas such as aircraft or helicopter 
landing, transfer of equipment or cargo at sea and off- 
shore installations. 

In the present study the estimation and prediction of 
heave, pitch, roll, sway, yaw motions of a DD-963 destroyer 
is considered, using Kalman filter techniques, for appli- 
cation in VTOL landing. 

The governing equations are obtained from hydrodynamic 
considerations in the form of linear differential equa- 
tions with frequency dependent coefficients. In addition 
non-minimum phase characteristics are obtained due to the 
spatial integration of the water wave forces. 

The resulting transfer matrix function is irrational 
and non-minimum phase. The conditions for a finite- 
dimensional approximation are considered and the impact of 
the various parameters is assessed. 

A detailed numerical application for a DD-963 destroyer 
is presented and simulation results of the estimations 
obtained from Kalman filters are discussed. The effect of 



2 


the various modeling parameters on the rms error is 
assessed and simplifying conslusions are drawn. 

The models developed are used to predict the motions 
a few seconds ahead. An upper bound for prediction time 
of about five seconds is established, with the exception 
of roll which can be predicted up to ten seconds ahead. 
The effect of noise and modeling errors on the rms pre- 
diction error is investigated in detail. 
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INTRODUCTION 

The present study started as part of the effort 
directed toward designing an efficient scheme for landing 
VTOL aircraft on destroyers in rough seas. A first study 

showed a significant effect of the ship model used 
on the thrust level required for safe landing. 

In a landing scheme therefore it would be desirable 
to have accurate ship models capable of providing a good 
real time estimation of the motions, velocities and 
accelerations of the landing area, resulting in safer oper- 
ations and with reduced thrust requirements. 

The modeling is quite complex and a substantial effort 
is required to reduce the governing equations to a finite 
dimensional system of reasonable order. 

The study contains a first chapter on the equations 
of motion as derived from hydrodynamics, their form and 
the physical mechanisms involved and the general form of 
the approximation. 

The second chapter describes the modeling of the sea, 
which proved to be a crucial part of the overall problem. 

The third chapter describes the derivation of the 
state-space equations for the DD-963 destroyer. 

In the fourth chapter the Kalman filter studies are 
presented and the influence of the various parameters is 


assessed. 
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In the fifth chapter the feasibility of predicting the 
ship motions a few seconds ahead in time is studied within 
the present formulation. 

Finally the appendices provide the characteristics of 
the destroyer, hydrodynamic information and some computer 
programs used. 
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OVERVIEVJ 


The real time estimation of the rigid body motions, 
velocities and accelerations of a vessel in rough seas requires 
accurate modeling of the wave exciting forces and the hydro- 
dynamic coefficients of the ship. 

The wave forces are obtained after an integration over 
the ship hull of the pressure forces, so that their evalua- 
tion requires a seakeeping program, while their magnitude and 
phase represent clearly an infinitely dimensional system with 
non-minimum phase characteristics. 

The complexity of the resulting equations is due pri- 
marily to the wave formation as the vessel moves, which is 
a mechanism of energy dissipation and additionally it introduces 
memory effects. 

The wave spectrum contains a rather narrow band of fre- 
quencies so that an efficient approximation of the ship charac- 
teristics can be achieved within this frequency band. 

A DD-963 destroyer was used as the basis for the present 
study. First the geometric and mass properties of the vessel 
were analysed by the M.I.T. Ocean Engineering Department Sea- 
keeping program and its hydrodynamic forces and coefficients 
were obtained. 

Subsequently a finite dimensional approximation was fitted 
in this data within the wave frequency range. Two groups of ship 
motions were distinguished, the heave-pitch and the roll-sway-yaw 
sets of motion, which up to the first order are uncoupled to each 
other. 

The parameters of the approximations are four: 
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The speed of the vessel, the direction of the waves, 
the significant wave height and the modal frequency of the 
wave spectrum. 

These models were used to estimate the ship motions, 
velocities, accelerations using noisy measurements of the 
motions. The Kalman filter designed for this purpose gives 
very good results when a relatively accurate estimate of 
the modal frequency of the spectrum is available. The 
modal frequency was found to be the most significant para- 
meter in the overall scheme since it influences the estimation 
error significantly and is the most difficult to estimate. 

The ship speed and the wave heading are important para- 
meters also, but can be estimated easily and accurately. 

The double peak spectrum, i.e. seas containing swell also, 
require separate treatment, because the low frequency peak is 
hard to estimate, while its influence is quite important. 

The predictability of ship motions has been investigated 
within the frame of the present study. First perfect state 
information is assumed and by propagating the prediction error 
covariance from zero initial value it has been established that 
within 25 % rms error over rms motion, the prediction time is 
about five seconds for all motions with the exception of roll 
which can be predicted up to ten seconds ahead. Simulations 
confirmed these results. 

The effect of noise and modeling errors is to reduce the 
prediction time. Omission of the non-minimum phase zeros has 

a particularly pronounced effect. 

In summary, the approximations described in the sequel 

provide a good model of the quite complex ship equations of 
motion within the wave frequency range. The derived models can 
be used for a real time estimation and prediction of the ship 
motions and other responses using Kalman filter techniques. 

Computer programs have been prepared that provide the 
required model matrices once the parameter of the problem has 
been specified. 



Chapter 1: EQUATIONS OF MOTION 


Definitions 

The rigid body motions of a ship in six degrees of freedom 
are shown in Figure 1.1; We define the plane to coincide with 

the symmetry plane of the ship, with the axis pointing vertically 
upwards when the vessel is at rest, and the y^ axis so as to ob- 
tain an orthogonal right-hand system, while the origin need not 
coincide with the center of gravity. The x^y^z^ system is an 
inertial system with x^y^ fixed on the xindisturbed sea surface, 
while the x y z system is moving with the steady speed of the 
vessel (i.e. it follows the vessel but it does not participate 
in its unsteady motion). Then the linear motions along the x^^, 

^1' ^1 surge, sway and heave respectively. In order to 

define the angular motions, we normally require Euler angles, rn 
the present case, though, we consider small motions so that the 
tensor of angular displacements can be replaced by a vector of 
small angular displacements, which are roll, pitch, yaw around 
the x^, y^, z^ axes respectively. 

The characteristics of a ship are its slender form, i.e. 

L/B>>1, L/T>>1, where L is the length, B the beam and T the draft. 
Also, the ship is symmetric about the xz plane and near symmetric 
about the yz plane. For this reason 


I 

I 


yz 

xy 


I 

zy 
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yx 


0 
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The value of 1 ^^ is typically small compared with 
The justification of using the linearity assumption is as 
follows: The excitation consists of wave induced forces, which 

include fluid inertia forces and hydrostatic forces. It is 
well established that the wave height to wave length ratio is 
small, since at a typical upper value of 1/7 the wave breakes 
and loses all its energy [15] (Figure 1.2). As a result, the 
major part of the wave force is a linear function of the wave 
elevation and can be obtained by a first order perturbation 
expansion of the nonlinear fluid eguation, using the wave height 
to length ratio as the perturbation parameter [15]. 

The wave spectrum, as will be shown later, has a frequency 
range between typically 0.2 and 2 rad/sec. Given the large mass 
of the vessel, the resulting motions, within this frequency 
range, are of the order of a few feet, or a few degrees, so 
that the equations of motion can be linearized. 

The only motion that requires attention is roll, because 
due to the slender form of the ship, the rolling motion may 
become large, in which case nonlinear damping becomes important. 

Simple Derivatibn 

We derive the equation of motion for a simple two dimensional 
object to demonstrate the overall procedure. 

Let us assume that we wish to derive the motion of a two 
dimensional cylinder subject to wave excitation, allowed to move 
in heave only (Figure 1.3). 
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The incoming wave of amplitude a^ and frequency will 
cause a force on the cylinder, and, therefore, heave motion. 

Due to the linearity of the problem, the following decomposition 
can be used, which simplifies the problem considerably. 

(a) Consider the sea calm and the ship forced to move 
sinusoidally with unit heave amplitude, and frequency and 
find the resulting force. 

(b) Consider the ship motionless and find the force on the 
cylinder due to the incoming waves and the diffraction effects 
(diffraction problem) . 

(c) In order to find the heave amplitude, within linear theory, 
we equate the force found in (a) times the (yet unknown) heave 
amplitude, with the force found in (b) . (Figure 1.4) 

The force in (b) can be decomposed further for modeling 
purposes, again due to linearity: One part is due to the un- 

disturbed incoming waves and the other part due to the diffracted 
waves. The first is called the Froude-Krylov force and the 
second the diffraction force . The total force is called the 
excitation force [15] . 

The force in (a) due to linearity can be also decomposed: 

The first part is simply the hydrostatic force. The second part' 
is the dissipative force, caused by the fact that the refraction 
waves carry energy from the ship to infinity. For this reason, 
we define a damping coefficient B so that the dissipative force 

I I 

will be -Bx where x is the heave velocity. The third part is 



16 


is an inertia force, caused by the fact that the heaying ship 
causes the fluid particles to move in an unsteady motion so that 

tl 

we define and "added” mass A and the inertia force becomes -Ax 

II 

with X the heave acceleration. If we denote the undisturbed in- 
coming wave elevation amidships as n(t) : 

n(t) = . ( 1 ) 

Where the real part of all complex quantities is meant, here 
and in the sequel. Then the excitation force will be 

F ^ ao • 

Where F^ is complex (to take into account the phase difference 
with respect to the wave elevation) , and the equation of motion becomes 

" ■ II ' I 

Mx = F - Ax - Bx - Cx (3) 

Where M the mass of the cylinder; the motion is also sinusoidal 
so with Xo complex; 

x(t) = XoS^“»^ • (4) 

A very important remark is that ^,A,B depend on the frequency 
of the incoming wave tOo • This can be easily understood by the 
fact that at -various frequencies the heaving cylinder will produce 
waves with different -wavelength. We rewrite, therefore, equation 
(3) as: 

-MXeWoS^^®^ = {Fo'(Wo)a^+ [A(w^)o)^ - im^ B(Wp) - C] Xo>e^“'»^. (3a); 

By dropping we can rewrite equation (3a) as: 


{-[M + A(a)o).l(/o + itOoB(Uo) + C} Xo = Fo(Uo) 


(3b) 
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The motion of a cylinder in water, therefore, results in an 
increase in the mass and damping term. Equation (3b) is used 
because of its similarity to a second order system, it is strictly 
valid, though, only for a monochromatic wave. 

Ultimately, we wish to obtain the response in a random sea, so 
equation (3b) must be extended for a random sea. This can be done 
by obtaining the inverse Fourier transform of (3a) , i.e. 

CO ■ ' 00 . ■ 

-i ^ - t) X (t) dT + 

+ C x(t) =_/ K^(t - T)n(t)dT 

VJhere the inverse Fourier transform of [M + A(oj)], 

iu)B(aj) and F„ (oj) respectively. The random undisturbed wave elevation 
is denoted by n(t). Equation (5) is not popular with hydrodynamicists , 
because the effort required to evaluate the kernels K^, K^, is by 
' far greater than that required to find the added mass, damping ana 
excitation force. For this reason, equation (5) is rewritten in a 
hybrid form as follows : 

-{M + A(w)] X(t) + B((jj) x(t) + C x(t) = F(w)n(t) (6) 

This is an integro-dif f erential equation (or differential 
equation with frequency dependent coefficients) , whose meaning is 
in the sense of equation (5) . 
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Strip Theory 

The evaluation of A(o)) , B((o), F(w) is not an easy task for 
complex geometries, such as the hull of a ship. The hydrodynamic 
particulars can be found in a later section, but we can give a 
simple description here of a technique used to simplify the 
derivations; [15], [17] 

The ship can be divided in many transverse strips as shown 
in Figure (!• 51 Due to its elongated form and for high frequencies, , 
each strip has small interactions with the other strips, except 
near the ends. Usually these end effects are small, so that 
instead of solving the overall three dimensional problem, we can 
solve many two dimensional problems (one for each strip) and 
sum up all the partial results. For the case of heave, for example, 
if A{w,x) B(oj,x) are the added mass and damping in heave of a strip 
at location x, then 

■ L/^ 

A(u) = / A(o),x) dx (7) 

-i/2 

i/z 

B(u) - / B(w,x) dx (8) 

-i/z 

The strip theory has larger errors at smaller frequencies. 

It so happens, though, that at small frequencies the hydrostatic 
forces are predominant, so that the motion error is quite small. 
Comparison with experiments has shown that for slender ship 



19 


configurations, the strip theory provides very good predictions 
[15] , [17] . 


Relation Between Added Mass And Damping 

The added mass and damping coefficients are not indepen- 
dent of each other, because their frequency dependence is caused 
by the saune refraction waves- If we define 

T(o)) = 0)2 [A(o)) - ^^] (9) 

0) 

Then ,T(a>) is an analytic function [16 ] • As a result, A(oj), 
B(w), which are real, are related by the Kramers -Kronig ■ relations, 
in order to describe a causal system. This fact will be used 
later to obtain a single approximation for T(o)) instead of two 
separate approximations for A(u)) , B(aj). 

Speed Effects 

As it can be seen in Figure 1.6 when the ship is heaving with 
a small angle 9 and at the same time is moving forward with speed 

O 

U, then a heave velocity results, which is x = U9- The effect of 
the forward speed, therefore, is to couple the various motions 
by speed dependent coefficients. As it can be found in Appendix 1, 
there are simplified expressions for the added mass, damping and 
exciting force v;ith a parametric dependence on the speed U. Then 
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QxprGSsions gxeatly facilitatG the evaluation of the ship motions. 

Frequency Of Encounter 

An additional effect of the ship speed is the change in the 
frequency of encounter. If the incident wave has a frequency oi 
and a wave number k, then the frequency of encounter is 

0 )^ = u + k U cos 4) (10) 

Where 4i is the angle between the x axis of the ship and the 
direction of wave propagation (Figure 1.7). In deep water, the 
dispersion relation for water waves is 

0 )^ = kg (11) 

so that we can rewrite (10) as 

2 

(ji = w +— U cos 4> (12) 

e g 

A very important consideration in the difference between 
frequency of encounter and wave frequency is the following. The 
ship motions due to linearity ;'7ill be of frequency so that the 
refraction waves are of frequency oi and the added mass and damping 
can be v;ritten as A(u ) , B(u) ). 

The amplitude of the exciting force though, consists of 
the Froude Krylov part which depends on u and the diffraction 
and speed dependent parts which depend on The time dependence 

is again w^’t, i.e. 
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F(t) = aoF(o), Wg/V) (13) 

with Eo the incident wave amplitude. 

This is a very crucial observation and can cause significant 
errors if not taken into account. 
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Equation of Motion 

Following the notation of Appendix 1, we write the equations 
of motion. It should be noted that, due to the slenderness of the 
ship, the surge motion is left out as a second order motion. This 
is in agreement with experiments [ VJithin linear theory and 

using the ship symmetry, the heave and pitch motions ar e not coupled 
with the group of sway, roll, yaw motions . This is not to imply 
that the motions are independent, because they are excited by the 
same wave, so there is a definite relation both in amplitude and 
phase . 

(1) Heave - Pitch Motions 


]) L + [ 1 

Oly A53A55 -V ] 


r C33 C35 -j 
‘■.C53 Css 


X = [ 
-V ^ 
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inhere A ^ j/ B ^ j,' 
efficient matrices; 


The frequency 
but is understood. 


...the added mass, damoing 
Fj the exciting forces; n 

= {X3, Xs)"^ 

X = {x2 , X4 / Xe } 


hydrostatic co- 
the wave elevation; 

( 16 ) 

( 17 ) 


u 

and velocity dependence is not written explic._ly. 


as described in the previous sections . 
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Heave - Pitch Approxiniation 

We start with the heave and pitch motions approximation. 

As it is obvious from equation (14) , it involves two stages: 

(a) Approximation of the exciting force 

(b) Approximation of the added mass and damping coefficients 

Data are provided by the hydrodynamic theory for both compon- 
ents and within the wave frequency range. 

A. Exciting Force Ap^proximation 

Figure 1.8 shows the exciting heave force on a box-like ship 
[25] . This information is important to demonstrate several zeros 
of the amplitude of the heaving force. Figure 1.9 shows the ampli- 
tude and phase of the exciting force on a destroyer, where, again, 
the same zeros appear, accompanied in the phase plot by jumps in 
the phase. 

The transfer function between the wave elevation and the heave 
force cannot be represented as a ratio of polynomials of finite 
degree as evidenced by Figure 1.9. Similar plots can be obtained 
for the pitch moment. Within the wave frequency range, though, only 
the first zero is important, while the remaining peaks are of minor 
significance. This is not true for other types of vehicles such as 
the semi— submersible , but for ships it is valid for both heave force 
and pitch moment, so it will be used to simplify considerably the 
modeling procedure. 

As it was mentioned before, the exciting force changes with 
frequency but its amplitude is determined on the basis of the 

frequency w. The following variables must be included in an 
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appropriate modeling of the exciting forces 



( 1 ) 

frequency w 



( 2 ) 

speed V 



(3) 

wave angle (}) 


F 3 (t, a„, <j). 

U) = 

F 3 (w, (j>)ao e « 

(18) 

Fs tt, Oo, <!>/ 

U) = 

{Fs (W, 4>) + f 3 (Wg,t)^ ® 

(19) 


where a, the wave amplitude, £3 the heave diffraction force. Equa- 
tions (18) , (19) show that the heave force does not depend on the 
ship speed, whereas the pitch moment does, in a linear fashion. 

In order to approximate F 5 (w,$) , f 3 (Wgf<t^ ) 

plots in Figure 1.9 as well as Figure 1.10, which show the approx- 
imate influence of the wave angle on the excitation force. 

In order to model the DD963 destroyer, the M.I.T. five degrees 
of freedom seakeeping program \21 ] was used to derive hydrodynamic 
results. The following model was derived to model shape of the 
heave force at V = 0 and (j) = 0 (no speed, head seas) 


Fa(s) 


tti n 



'( 20 ) 


Where J = 0.707,ai a constant to be determined from hydro- 
dynamic data, n the wave elevation and the corner frequency. 
Remembering the analysis above concerning the dependence of the 
force on u and Figure 1.10, we can derive „ 
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2 -n q 


Lcos(j) + B 


2 7T 


LCOS4> + B 


U COS<j) 


( 21 ) 


where L is the ship length, B the beam. 


Before we establish a relation similar to (20) above, we have to 
discuss Figure 11, where it is shown that for long waves, the heave 
force and the pitching movement are 90° out of phase . This means 
that the transfer function between heave and pitch is a non-minimum 
phase one, because the amplitude is constant, while the phase is 
90°. We choose to attribute the non-minimum phase to pitch. Also, 
the pitch angle tends to the wave slope for large v/avelengths , so 
the pitching moment can be written as 


F 5 - az 


1 - s/t»), 


cos4> 


7 — T 


1 + S/Wo [1 + 2J-f- + J J 

8i B. 


r\ 


( 22 ) 


V7here aj a constant to be determined, is the same (for 
simplicity) as in equation (21) and is an artificial frequency 
to model the non-minimum phase- It will be chosen to be equal to 
the wave spectrum modal frequency, so we defer the discussion until 
the corresponding section. 


B. Added Mass and Damping 

By using equation (9) , we can rewrite the equations of motion 


T 3 3 

s^ 

+ MS^ + 

C 3 3 

T35 

+ 

UST33 + 

C3 5 

Ts 3 

sV 

- UST33 

+ C3 5 

Is^ 

+ 

T55 s’- 

U^T3 3 + C55 


[pin (23) 

F 5 
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Here we construct a simplified model where 


Tij= Aij- 


withAij,Bij to be evaluated from the hydrodynamic data, 
tion (23) can be written, after we define 


Yi = X3 


ya = X3 


Then eaua- 


in the form: 


Ys = X5 


Y4 = X 5 


y = {yi /Y 2 /y3/Y4)^ 


A3 3 A3 5 1 

I 

A53 As 5 


C33 B33 C35 B35 
C53 B53 C55 Bss 


(23a) 


where 


C35 = C35 + U A33 
C53 “ C35 “ U A33 


C^s = C 55 - ^ A; 
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Sway-Roll-Yaw Approximation 

Next we approximate the sway-roll-yaw group of 
motions, which is uncoupled to first order from the 
heave-pitch group of motiqns . Again a two stage app- 
roximation is required, i.e.; 

(a) Approximation of the exciting force 

(b) Approximation of the added mass and 
damping coefficients 

Data are obtained by using the sea-keeping program. 

A. Exciting force approximation 

The same infinite-dimensional form is obtained for 
the exciting force as seen in Figure 1.12 (17abc) for all 
three motions, as in the case of pitch and heave. Again, 
within the wave frequency range, a finite dimensional 
approximation can be achieved, and of reasonably small 
order. 

The important fact that the exciting force depends 
on the wave frequency rather than the frequency of en- 
counter, is used, while the following three quantities 
define the exciting force amplitude and phase 

(1) Wave frequency 

(2) Speed V 

(3) Wave angle cj) 

In Appendix 1 the strip theory approximation of the 
sway, roll, yaw forces can be found. Using the M.I.T. five 
degrees of freedom seakeeping program the following finite 
dimensional approximation was found in case of V = 0 , (|) = 
90® for sway, roll and yaw 
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F 2 (S) 


A2 


2 

S 



2 

+ 2J2 



1 


(27) 


F 4 (S) 


S 

W 4 


2 + 2 J 4 — + 1 

CO4 


(28) 


Fe (S) = 


s 

W6 


2 S 

+ 2Je ^ + 1 

We 


(29) 


where 


w = 
2 

0.65 

J = 
2 

.5 

(i) 4 = 

0.85 

J 4 = 

.3 

Wg = 

0.85 

Je = 

.3 (30) 

are 

obtained 

from hydrodynamic data. 


We redefine the value of W 2 f W 4 ,We such that 
it will be valid for angles (j) other than 90°, and 


(31) 


speed other than 0 : ^ 

w' . = (o) . + (jo? — cos 6 ) sin<j) 

3 D 3 9 

where j = 2,4,6 and is given above. 

It should be noted that the sway, roll and yaw forces 
are proportional to the wave slope, i.e. 90 out of phase 
with respect to the wave amplitude. This means that they 
belong to the same group with pitch, and the same non- 
minimum phase transfer function; 


s - Wq 
s + Wo 

must be used for all three of them, when the total 
system (all 5 degrees of motion) is considered. 
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B. Added Mass and Damping 

The amplitude of the transfer function between the wave 
elevation and the rolling motion has a very narrow peak so that 
the coefficients can be approximated as constant [17], Using 
Appendix 1: 


N‘4 

= A° 
44 

®44 

®44 

*42 

= A® 
^24 

^42 

= ®24 

*46 


V 

"2^24 


®46 


- V A® 

^ ^24 



using the value of w at the roll peak. It should be noted 
that roll involves a significant nonlinear (viscous) damping, 
which is approximated by introducing an additional "equivalent 

damping coefficient [ 3 1 
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Eiirdlarlv, we calculate the sway, yaw coefficients at 
the same frequency: 


'22 


= A 


22 


22 


= b; 


A 


26 


= A 


26 


;^22 


B 


26 


B 


26 


- V A 


22 


62 


B 


62 


A 


66 


B 


66 


A 0 _ _ B ® 

22 


B.“ + V A° 


26 


22 


66 ^2 22 


oO + Y. B^* 
®66 „.'"22 


(33) 


SMA..+Mij} + S{B..) + {C..) 


A 


X 



i = 2,4,6 
j = 2,4,6 


(34) 


where . = 0 except for which is the roll hydrostatic 

constant, i.e. = A* (GM) with A the ship displacement and (GM) 

the metacentric height. 

Due to the special form of the matrix C, a zero-pole 
cancellation results from a direct state space representation of 
the equation above. After some manipulations, the following 
representation can be obtained which avoids zero-pole cancellation 
problems : 
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X = T X + U F 


where X = {Xg / , X^ , Xg} 


F = {/F^, F^ , F^ , /F,, F J 


where /F indicates the time integral of F and T = (t-}and 


U = {u. . } , with; 
r 




. X 


12 


r r 32 r r 

U\ = - -p , t42=;r- ' t43= :r^^P - P / t44 = T^^P - p 

r22 21 31 ^22 ^22 22 32 ^22 23 33 


■»- =— P -I- — P t = — P t — P t — P 

^21 21^11 23S1 22 2 1 12 2 3 42 22 


^2 3 ^ 21 ^ 13 ” ^ 23 ^ 43 “ ^2 2^4 4 


t = “P t -P t 

H4 2 1^14 23S4 


tjj^ = 0 except for ~ ^ 


“ij 


’=12 ’=21 


= 0 except for ^ 


22 


U. 


14 


. .r , r 

12 23 

“^13" ~F“ ' 
22 


U = r 


,r . . r 

32 2 1 


■*1 3 1 

22 


r r 

u = r - , 


44 331 ^ 


22 


U 21 = -P u -P u 

2111 2341 


U =r,U =r,U =-PU-PU 

22 21 23 22 24 21 14 23 44 


^25 ^2 3 


where R = = [A+M] ^ 


(35) 


(36) 


(3.7) 




P = {P^j} = R B 


(38) 
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Figures 1.13 through 1.15 show the actual force and moments 
versus frequency and the achieved finite dimensional approximation 
for zero speed and beam seas (({)=90°). It should be noted again 
that the approximation is not as good outside the wave frequency 
range. 

Figures 1.16 through 1.18 show the same quantities for speed 
U= 15.5 ft/sec and 45° angle of incidence. From these figures it 
can be seen that above 1 rad/sec the approximation is poor, none- 
theless no significant wave energy is contained in that range, so 
the approximation is acceptable. 

Figures 1.19 through 1.21 show the overall transfer function 
between the corresponding motion and the sea elevation. 
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5^A »i K/> Vr 


Figure 1.2 
Wave Profile 



Figure 1 . 3 


Hydrodynamic Problem 
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c^Qitvcic.r 
p<i SI Vu»n 


Figure 1.4a 
Radiation Problem 



Diffraction Problem 



37 


Figure 1.5 

Strip Subdivision of Ship Hull 





Figure 1.6 

Effects of Ship Speed 
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^* 180 * 






Fn= V/ygT 


Figure 1.9 

Effect of Ship Speed on Heave Force and Pitch Moment 

at Head Seas 



Figure L.IO 

The Equivalent Wavelength for Heave and Pitch 
is X/cos<t) 



Figure 1.11 

Phase Difference Between Heave and Pitch 

\ 
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Figure 1.12 

Wave Force and Moments in Sway, Roll, and Yaw at 
60° Angle and Zero Forward Speed 
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Figure 1.13 

Sway Force Versus Frequency for the DD963 Destroyer and Its 
Finite Dimensional Approximation (dotted line) for Speed 
U = 15.5 ft/sec and 45® Angle of Incidence 
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Figure 

1.14 



DD963. 


Same conditions as in 1.13 




Rmp 1 1 tude (dB) 


45 


SWRY FORCE 



Omega (dB) 

Figure 1.16 

Sway Force Versus Frequency for the DD963 Destroyer and Its 
Finite Dimensional Approximation (dotted line) , for Speed 
U = 15.5 ft/sec and 45° Angle of Incidence 
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Figure 1.17 


Roll Moment for DD963. Same conditions as in 1.16. 
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Figure 1.18 

Yaw Moment for DD963. Same conditions as in 1.16, 
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ROLL TRRNSFER FUNCTION 



Ome (dB) 

Figure 1.20 

Roll Transfer Function for DD963. Same conditions as in 1.19. 
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Chapter 2 : SEA MODELING 


The sea waves are generated by the wind, except for very few 
cases. The process of wave generation is of importance in model- 
ing, so we will outline, briefly, a simple theoretical model: 

When the wind starts blowing over a calm water surface, it 
contains gust components of high frequency, which cause wavelets 
on the surface. This is due to the inherent instability of the 
wave air interface. As soon as the surface becomes rough, a 
significant drag force develops between air and v/ater, which 
becomes zero only if the average wind speed (which causes the 
major part of the drag) equals the phase wave velocity. As a 
result, the steady-state condition of the sea develops slowly 
by creating waves whose phase velocity is close to the v;ind speed. 
Since the process starts with high frequencies, we conclude that 
a young storm will contain a peak at high frequency. We usually 
distinguish between a developing storm and a fully developed 
storm. 

As soon as the wind stops blowing, then the water viscosity 
dissipates the high frequency waves so that the so called swell 
(decaying seas) forms, which consists of long waves (low frequency 
content) , which travel away from the storm that originates them. 
For this reason, swell can be found together with another local 
storm (Figure 2.1). 
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A storin usually contains one peak (except if swell is present 
when it contains two peaks) and the peak frequency is called 
the mod al frequency (Figure 2.2). AlsO/ the intensity of the storm 
is required, which can be described in a number of waves; Beaufort 
Scale, Sea State, Wind Average Velocity, Significant Wave Height. 

The best is the significant wave height H defined as the statistical 
average of the 1/3 highest waveheight. For a narrow band spectrum 
of area Mo 


H = 4 / Mo 

From our discussion on sea storin generation, we conclude that 
it is important to model a storm by both H (intensity) and 
(duration of storm) . For this reason, the Bretschneider Spectrum 

will be used defined as: 


S(m) = 


1.25 


Oi 






m 


0 ) 


exo { - 1.25 (■ 


m 


u 


■r } 


(40) 


The spectrum was developed by Bretschneider for the North 
Atlantic, for unidirectional seas, with unlimited fetch, infinite 
depth and no swell. It was developed to satisfy asymptotic 
theoretical predictions and to fit North Atlantic data. It was 
found to fit reasonably v/ell in any sea location. Also, by 
combining two such spectra, we can model the swell as well. Its 
main limitations are unidirectionality and unlimited fetch. 
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It was felt, however, that it could provide an adequate descrip- 
tion for the present application for open sea. 


As it has already been mentioned, the forward speed of the 
vessel causes a shift in the wave frequency to the frequency of 
encounter. The dpectrura, now, can be defined for ship coordinates 
as follows 


S(^^) 


S(M) 

doj /do 
e 


(1) = f (Wg) 


(41) 


where 


w = f (O)^) 


-1 + /T 


+ 4o). 


Ucosi 


2 


— cosA 

g 


(42) 


A rational approximation was found to (29) subject tp (30) in 
the following form 


S»(“ ) = 

a e 


1.25 

4a)^ 

m 


B(a) 


(o) /o) ) 
e' o 


4 


1 + ( 



(43) 


where S (w ) the approximate spectrum 

Si 0 


a 



63 COS(b 

in 


( 44 ) 


B(a), y{a) functions given in table 2.1 
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Now a transfer function can be deriven from (43) such as to 
provide an output with the spectrum in (40) when driven by white 
noise. It is easy to see that 


H^(s) = 

a 


(s/w^)^ 



(45) 


where 



1.25 


4u 

m 




B(a) 


(46) 


Wo = Y(a) 


(47) 


J = 0.707 

A plot of the spectrum for various wind speeds is given 

in Figure 2.3, while the spectrum and its rational approximation 

1/3 

is plotted in Figure 2.4 for fully developed seas and H = 3 m. 


Important Remark 

It is customary to define the power spectral density as the 
Fourier transform of the autocorrelation R(t) 

S(w) = f R(t)e dr (48) 

— 00 

In wave theory, the spectrum is defined one sided (for positive 
frequencies only) as follows 


S^tw) = 


/ X -IWT 
/ R(T)e dx 


CO 


> 0 


(49) 



55 


For this reason, the relation between the spectrum S(u) as 
required for the present application and the Bretschneider spectrum 
Sg(o)) is: 


S(o}) 


r 

ir Sg(io) 

TT S„ (-u) 


0 ) > 0 

(I) < 0 


(50) 


Therefore, the intensity of the white noise required for 
driving the transfer function (33) is r (or equivalently we can 
multiply the transfer function by /W”) . 
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TABLE 2 . 1 

Sea Spectrum Coefficients 


a 

. 00 
. 10 
.20 
. 30 
. 40 
. 50 
. 60 
.70 
. 80 
.90 
1 . 00 
1 . 10 
1 . 20 
1 . 30 
1 . 40 
1 . 50 
1 .60 
1 . 70 
1 .80 
1.90 
2.00 


2 . 7306 
2.8218 


Y ( a ) 

3 (ot) 

. 9538 

1 . 8861 

1 . 0902 

1.6110 

1 . 1809 

1 . 3827 

1.2717 

1.2116 

1 . 3626 

1 . 0785 

1 . 4539 

.9718 

1 . 5448 

. 8845 

1 . 6360 

.8116 

1.7272 

. 7498 

1.8182 


1 . 9095 

. 6509 

2. 0008 

. 6106 

2.0918 

. 5750 

2.1833 

. 5434 

2.2744 

. 5150 




: S 

2.5481 


2.6395 



. 4085 












eii seas 



„ OJ 

cui/iec) 



PECfRflL DENSITY FUNCnON 
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F R EQU E N C Y (rad /' s e -s ) 


Figure 2.3 Bretschneider spectrum at various 
wind speeds (fully developed seas) 
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Chapter 3:; DERIVATION OF THE STATE - SPACE EQUATION 

We proceed to derive a state-space form of the 
equations of motion. Starting with the sea, we can construct the 
following representation {three cascaded second order systems) 

1 0 0 0 0 

-2Jojo 0 Wo 0 0 

0 0 1 0 0 

-s 

0 -a?o -2JWo 0 w^o 

0 0 0 0 1 

000 -Wo “2JWo 

n=[ -y^OOOOO] Xg 




or 


o 



""s ^s 


+ B w 
s 


(52) 


Hea ve-Pitch Model 

The following model is derived for the force (some algebra 
involved to reduce the dimension of the state) • 


was 
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where 





Di 

D2 


A 3 3 

A 3 5 

” 1 








D3 

D„ 

• 

A 53 

A 5 5 





B^ 

B® 

B^ 



Di 

D2 


C 3 3 

B 3 3 

C 35 

B 3 5 



B® 

B’ 

B® 



D3 

D4 


C 53 

B 5 3 

C 5 5 

B 5 5 


( 56 ) 


(57) 


and we can write in short 

>r 

Fs_ 

(58) 

= C X.m 
m — 


X3 

X5 


X„ = X 
-m m -m 


B 


m 


Now the total model can be constructed as follows; 


xt = \ ^t ^ ®t ^ 


(59) 


X3 
X s 


= 


where 


r 


^t = 


A 


A B ! 
m m f t 


B.C 0 
f s 


u 


^ = 


B 


(60abc) 




o 


c o 

m 
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Data were obtained for the DD-963 destroyer from the M.I.T. 
Seakeeping Program and are given in Appendix 2. Appendix 4 lists 
a computer program that produces the A^, matrices once the 

ship speed, wave heading, significant wave height and modal 
frequency were specified. The output can be used directly as 
input to the LIDS control and filter design package. 

Table 3.1 provides the numerical values of matrix A for 
speed 20 knots, angle 0° (head seas) , significant wave height 
10 ft. and modal frequency 0.72 rad/sec (sea state 5). 
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Sway-Roll-Yaw Model 

The sway, roll, yaw exciting forces and moments are essentially 
driven by the slope of the sea elevation, which for regular waves 
equals the wavenumber times the amplitude, or for deep water we 

can write 

Slope = 


i.e. in the time domain: 

Slope = 


d 
dt^ 


(61) 


As outlined in Chapter 1, the fact that only roll has a, 
spring constant causes zero pole cancellation problems, whxch can 
be avoided by introducing the matrices T and U described by 
equations (36) through (38). The original equation is in the form 

(M+A) x^ + Bx^ + Cx^ “ Hi 


where x^ 


is of dimension 3: 


^1 = 


' 

sway 


■xr 

X4 

Xs 


'F 2 ' 

F 4 

roll 

yaw 


II 

Fg 


(63) 


M is the mass matrix, A and B the added mass and damping 
matrices respectively, C is the hydrostatic matrix and F the vector 
of exciting force and moments. Then 

S^Xq = - [M+A]"^ sx^ - [M+A]"^Cx^ + IM+A]"^F^ 


(64) 
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By letting 

T ' 

X = [X2 , X4 , X4 , Xe ] 


f'^ = [^F2, F2, F4, ^F4, Fe ] (65ab) 


and using the T and U matrices of equations (36), (38) we obtain 
a state space description of (64) without zero-pole cancellation, 
in the form 

X = T X + U F (66) 


The state space representation of equations (27) through (29) 
is in the form 


0 

1 


0 

1 


L : 


X . + 


n 

(67) 

-2 ^ . u) . 

3 3, 

-3 

_1 



where j = 2,4,6, while 


^F.' 

3 


A 

3 3 

0 

F . 

II 

0 

A .03^. 

L 3 _ 



3 3 J 


( 68 ) 


For this reason we build a force matrix A^: 


^F 


A2 

0 

0 


0 

A 4 

0 


0 

0 

As 


(69) 


and a matrix driving the force dynamics with n , i.e. using the 

r 

sea model, which is exactly the same as in the case of heave, 
pitch (6 states) , so 
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Bp(i,j) = 0 

except 

B„(2,2) = -A2w| 

r 

B„(4,2) = -A 4 W 4 

r 

Bp (6,2) = -As We 


(70abcd) 


Then using the same sea model described in equations (52ab) we 


obtain the overall model as 


I 



(16x16) (16x1) 


X 2 I (3x16) 

X4 = x^ 

X6 


(71ab) 


where 



0 0 

Ap 0 
U T 


B^ = [ 0] 


10 0 0 
0 0 10 
0 0 0 1 


(72abc) 


Data for all quantities involved are given in Appendix 2, while 

Appendix 4 lists the computer programs that can produce the matrices 

A , B , C once the ship speed, wave heading, significant wave 
t t t 

height and modal frequency are specified. 
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Table 3.1 provides the matrix for speed 15.5 knots, 
heading 45°, significant wave height 10 ft. and modal frequency 
0.72 rad/sec (sea state 5). 



TABLE 3.1 


Matrix A for heave and pitch (15 x 15) 


1 . 

491E+08 -1.726E+00 

a. 

a. 

0. 
a. 
a. 
a. 
a. 
a. 
a. 
a. 
a. 

ae a. 
a. 


a. a. 

0. 1.491E+00 

0 . 1 . 

-i.49ie+B0 -'i.72se+aa 

a. a* 


a. 

a. 


a. 

a. 

a. 

a. 

0 . 

0. 

0. 

0 . 

a. 


a. 

a. 

a. 

a. 

a. 

a. 

a. 


0. 1.491E-i-a0 

a. 1. 

--I.491E+0a ~i.72BE+aa 

a. a. 


a. 

a. 

a. 

a. 


a. 

a. 

a. 

0» 

0. 

1 . 


a. 

a. 

a. 

a. 

a. 

a. 

a. 


a. 

a. 

a. 

a. 

a. 

a. 

a. 

a. 


~1.2l3E+00 -5.350E~01 


0. 

8.818E-04 

a. 

a. 

a. 

a. 

0 . 


a. 

3.106e-04 

a. 

a. 

a. 

a. 

a. 


i a. a. a. 0- 

a. a. 

a. a. 

a. 0- 

a. a* 

a. a. 

a. a. 

-9.452E+08 -2.190E+01 7.850E:“01 0. 

a» K a* 0* 

-1.1B3E+00 “4.347E-01 “1.529E-03 0. 

!. 


a. 

a. 

a. 

a. 

a. 

a. 

a. 


a. 

a. 


a. 

a. 

a. 


a. 


a. 

a. 

a. 

a. 

a. 


a. 

-B.B35E-01 -l.l52E+aa 

a. a. 


a. 

a. 

a. 

0 . 

a. 

a. 

a. 

a. 

a. 

a. 

a. 

S.835E-01 

0. 


a. 

a. 

a. 

a. 

a. 

0 . 

a. 

a. 

a. 

a. 

a. 

a. 

I. 


a. 

a. 


a. 

a. 

a. 

-S.027E-01 

a. 

1.60SE-a2 

a. 

a. 

a. 


a. 

9.72BE-01 


0. -B.S35E-01 -1-.152E+00 0. 

1.249E+00 -9.728E-0: 4.402E-01 ~1.050E+80 


Matrix A for sway, roll, yaw (15 x 16) 


GO 


1 . a. 

138E+00 -1.507E+00 0. 

0. 

a. 


a. 

1.138E+00 
0. i. 

~1.13SE+00 “1.507E+00 

a. a. 


a. 0 . 

a. a. 

5.722E+01 0, 

a. 0 . 

8.8S9E+02 0. 

0. 

6 


.461E4>03 


0. 

0. 

a. 

a* 

a. 

a. 

a. 

a. 

a. 

a. 


a. 0 . 

a. a. 

a. a. 

a. 1.13BE4-O0 

a. I. 

-1.13BE+00 -1.507E+00 


0. 

0. 

a. 

a. 

a. 

a. 

a. 

a. 


a. 

a. 

a. 

a. 

a. 

a. 

a. 


a. a. 

a. a. 

a. 0 . 

0 . 0 . 

0. a. 

a. a. 

0 . 1 . 

'2.6I0E--01 -5.202E-01 

0 . 0 . 


0. 

a. 

a* 

2.423E-03 

3.2IBE-08 

0. 

-4,0S7E-06 


0. 

0. 


0. 

1.375E-08 


a. 
a. 
a. 

0 . 

0. 

0. 

0 . 
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0 . a. 

0 . 0 . 

0. a. 


a. 

0 . 

a. 

a. 

a. 

a. 

a. 

1. 


a. 

a. 

a. 

a. 

a. 

a. 

0 . 

a. 

a. 

a. 

a. 


0 . 

0. 

0. 

a. 

a. 

a. 

1. 


0 . 

0. 


-8.111E-0I ^4.458E-01 
-4 . 469E-08 0. 

3,142E-06 -I.999E-08 -l.e04E^07 
8- 0. 0. 

0. J.299E-07 0. 


a. 

a. 

a. 

0 . 

0 . 

0 . 

0. 

0. 

0. 

0. 

a. 

a. 

-0.343E-03 

4.343E-05 

0. 

-4.6S1E-B4 


0. 

0. 

a. 

a. 


a. a. 

a. a. 

a. a. 

a. 0. 

a. a. 

2.740E-01 1. 

1.98IE-02 -2. 
1 . 0 . 
2.319E-02 4. 


a. 

0 . 

0. 

0 . 

0. 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

897E-01 8.041E4^00 

344E-01 1.2e6E-02 

0 . 

092E-04 -3.242E-02 
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Chapter 4 : KALMAN FILTER AND SIMULATION 

The heave-pitch approximation resulted in a 15 state system 
and the sway-roll-yaw approximation in a 16 state system. Given 
that 6 states describe the sea, the total system required for 
5 degree of freedom motion studies would contain 25 states. If the 
sea spectriam contains two peaks then a 31 state model is required. 

The heave-pitch group is not coupled with the sway-roll-yaw 
group so that the study of each group can be independent. This 
is not to indicate that in a total design the two groups must 
remain independent, since they are excited by the same sea. 

Heave-Pitch Motions 

It is assumed that the heave and pitch motions are measured. 

' The gyroscopes can provide accurate measurements of angles , 

up to about 1/10 degree. The noise therefore is due to structural 

vibrations, which in the longitudinal direction can be significant 

due to the beam-like response of the vessel. As a result the 

measurement noise was estimated based on data from ship vibrations. 

The same applies to the heave measurement noise. 

A Kalman filter was designed for speed V = 21^^/sec and waves 
coming at 0“ (head seas) with significant wave height H = 10 ft. 
and modal frequency = 0.73 rad/sec (sea state 5). The measure- 
ment noise intensity matrix was selected from ship vibration 
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data to be 


V 


0.75 

0 


0 

0.0003 


The model poles are shown in table 4.1, while the filter poles are 
within a radius of 1.3 rad/sec as seen in table 4.2. Typical 
simulation results are shown in Figures 4.1 and 4.2. In these 
figures exact knowledge is assumed for the significant wave height 
and modal frequency. The accuracy of the filter is very good both 
for heave and pitch. 

Subsequently, the same filter was used combined with a ship 
and sea model different than the nominal one, to investigate the 
sensitivity to the following parameters: 

The influence of the significant wave height is very small 
when the modal frequency is accurately known. On the contrary, 
the influence of the modal frequency is quite critical; particularly 
for pitch (see Figures 4.6 and 4.7). The same conclusion is 
reached when a double peak sea spectrum is used [22], [23]. 

The effect of the forward speed and wave direction was found 
to be unimportant particularly for heave, while for small changes 
in wave angle (+ 15°) the pitch prediction error was not affected 


significantly [22] . 
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Sway-Roll-Yaw Motions 

As in the case of heave and pitch, the measurement noise 
consists primarily of structural vibrations rather than instrument 
noise. For roll such vibrations are quite small for a destroyer 
vessel and similarly for sway and yaw the vibrations are smaller 

than in the case of heave and pitch. 

The noise intensity used was nonetheless similar to the 
heave, pitch noise, so as to bound the filter eigenvalues below 
2 rad/sec, which is the typical wave bandwidth. 

A specific example has been worked out for a forward ship 
speed of 15.5 ft/sec and waves at 45® and sea state 5 (significant 
wave height of 10 ft. and modal freqeuncy of 0.72 rad/sec). The 
measurement noise intensity matrix was 

diag {0.1 ft^ , 2 * 10 "“ (rad)2 , 2*10"“ (rad )2} 

The simulation shows very good estimation as seen in figures 4.3, 
4.4, and 4.5. Yaw is very small and the measurement noise is large 
relative to the yaw motion, nonetheless the yaw estimation based 

primarily on the roll, sway measurements is very good. 

Table 4.2 presents the results of a sensitivity study oj. i_he 

influence of the various parameters involved. The most critical 
parameter is again the modal frequency. The ship speed and the wave 
direction are not critical for the estimation error. This is a 
very important conclusion as far as the wave direction is concerned, 
because in reality, seas are directional and very difficult to 
measure, or even model appropriately. 
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The influence of systematic measurement errors was studied 
by using a calibration factor. This factor is defined to be the 
ratio of the measurement fed to the filter over the actual 
measurement, thus introducing a systematic error. If C is the 
calibration factor, then the systematic error as a percentage of 
the actual measurement is 100* (1-C). In the case of a 10% error, 
the most significant change was found in the case of the roll 
motion. In the case of a calibration factor 0 (indicating a dis~ 
connected measurement) significant errors resulted, especially 
for roll in the case of disconnected roll measurements (Table 4.3). 

Table 4.4 presents the poles of the model used, while Table 
4.5 shows the poles of the Kalman Filter derived for the nominal 
condition as described above. Figures 4.8 through 4.10 are sim- 
ulation results and show the significant effect of the modal fre- 
quency on the estimation error. Finally, Figure 4.11 shows the 
simulation of the sway motion estimation when the roll measurement 
is disconnected. 
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TABLE 4.1: Poles of the heave, pitch model 





74 


Ship Speed: U = 21 ft/sec. Heading angle: 0° 



TABLE 4.2 Poles of the heave, pitch Kalman Filter 






TABLE 4 . 3 


-T5 


Sensitivity of the RMS error of sway, roll, yaw motion, 
to changes in the parameters of the ship-sea model. 

C (sway) indicates a calibration coefficient in the 
sway measurement (and similarly for the other motions) , 
to handle systematic errors in the measurements. 


Parameter Changed 

Error Sway (ft) 

Error Roll (deg) 

Error Yaw (deg) 

Basic Case 

0.241 

0.56 

0.0776 

U=20 ft/sec 

0.245 

0.568 

0.0963 

Wj^=0.52 rad/sec 

0.314 

0.91 

0.0858 

4>=60° 

0.296 

0.624 

0.112 

C (sway) =0 . 9 

0.255 

0.586 

0.081 

C(roll)=0.9 

0.247 ' 

0.708 

0.0808 

C (yaw) =0 . 9 

0.242 

0.56 

0.0777 

C (sway) =0 . 0 

0.518 

1.21 

0.1408 

C(roll)=0.0 

0.376 

4.08 

0.158 

C (yaw) =0 . 0 

0.242 

0.563 

0.0785 

RMS values of 
the motions 
(nominal case) 

0.60 

4.56 

0.227 

Measurement 
noise intensity 

(.316 ft)^ 

(.81°)^ 

(.81°)^ 










TABLE 4.5: Poles of the sway, roll, yaw Kalman filter 
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Figure 4.1 

Results of Heave Simulation and Its Kalam Filter Estimate 
(dotted line) , Using Accurate Model at U=21 ft/sec and 0° 
Angle of Incidence and in sea state 5. 
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Figure 4.2 

Results of Pitch Simulation and Its Kalman Filter Estimate. 
Same conditions as in 4.1. 
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Figure 4.3 

Results of Sway Simulation and Its Kalman Filter Estimate 
(dotted line). Using Accurate Model at U=15.5 ft/sec and 45 
Angle of Incidence, and in sea state 5. 
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Figure 4.4 

Resulis of Roll Simulation and. Its Kalman Filter Estimate. 

Same conditions as in 4.3. 
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Figure 4.5 


Results of Yaw Simulation and Its Kalman Filter Estimate. 

Same condition as in 4.3. 
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Figure 4.6 

Results of Heave Simulation and Its Kalman Filter Estimate 
(dotted line) . The actual wave spectrum model frequency is 
0.52 rad/sec, while the value used in the Kalman Filter is 
0.72 rad/sec. All other parameters as in 4.1. 
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Figure 4.7 

Results of Pitch Simulation and Its Kalman Filter Estimate. 


Same condition as in 4.6. 
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Figure 4.8 

Results of Sway Simulation and Its Kalman Filter Estimate 

(dotted line). Actual =0.52 rad/sec, while in Kalman Filter 

m 

Wj^=0.72 rad/sec. All other parameters as in 4.3. 
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Figure 4.10 

Results of Yaw Simulation and Its Kalman Filter Estimate. 


Same conditions as in 4.8. 
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Figure 4.11 

Results of Yaw Simulation and Its Kalman Filter Estimate 
(dotted line) , using noisy measurements (light line) when the 
roll measurement is disconnected. Same other conditions as 
in 4.3. 
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Figure 4.12 

Results of Sway Simulation and Its Kalman Filter Estimate 
(dotted line), using measurements (light line) when the actual 
angle of incedence is 60° and the value used in the Kalman 
Filter is 45°. Same other conditions as in 4.3. 
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Chapter 5 : SHIP MOTION PREDICTION 

It is of interest to use the models developed in the previous 
sections to forecast the behavior of the vessel within a few seconds. 
The feasibility to predict the motions could assist significantly 
the pilot in committing the aircraft to landing under favorable 
conditions. 

An automatic landing does not require within LOG theory such 
information since the predictable part of the motions is included 
in the state and therefore used directly. Nonethless the prediction 
is of primary importance for pilot landing or semi-automatic landing. 

Similarly for offshore operations a display of a prediction 
of the most critical vessel motions could reduce the operation risk 
significantly. The operator could choose a time window of minimal 
motion or acceleration and then transfer cargo or personnel. 

The subject has been considered in the literature [6] ^ [HI r 
[24] using both frequency and time domain techniques. 

Theoretical Background 

The first to treat the subject of developing a predicting 
filter was Wiener [29]. If a random process has power spectrum S(w) 
a spectral factorization is first required, i.e. 

A(o)) = ip (w) i|;* (w) (73) 

where * denotes complex conjugation and ijj (s) is an analytic function 
of s with the exception of a finite number of poles in the left 
half plane. Then the transfer function of the optimal predictor 
K(o 3 ), in the sense of minimizing the expected value of the error, is 
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given by the expression; 


K(o)) 



, , -i(jJt 

Y(t+a)e 


dt 


(74) 


where Y (t) is the inverse Fourier transform of (oj) and a is the 
prediction time. The importance of this result is to provide a 
number of intuitive results, such as the fact that a narrow band 
process is predictable while at the extreme a wide band process is 
unpredictable. The disadvantage of this approach is that it may 
require differentiators in its implementation, depending on the 
form of the spectrum [29] . 

The alternative is to use state space models where no such 
problems appear. In fact the predictive filter has a very simple 
form. If the system has a state space description 


X = Ax + W 
Y = Cx + W 2 


(75ab) 


where and W 2 are white noise signals , it is not hard to see that 
if the state x is perfectly known at t then the predictable part of 
x(t+T) , denoted by x^, is 


X (t+T) = e'^^ x(t) 
^ — 


(76) 


If the state is not available the Kalman filter estimate is 
used instead [5], [8]: 


At 


X (t) 


X (t+T) = e' 


(77) 
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i.e. the propagation of the equation 


^ (t+T) = A X (t+T) 

✓s 

X (t) = X (t) 


(78ab) 


For a stable matrix AjXp(t+T) -> 0 as t ->-00, reflecting the fact 
that the influence of the driving white noise completely alters the 
state of the system, once the homogeneous solution has died out. 

The covariance of the error 


e = Xp - 

denoted as is governed by the equation 

O 


P (t) = A P (t) + Pp(T) ^ 


(79) 


When the state is perfectly known the initial condition is: 

P (o) = 0 (80a) 

It 

While in the case of using the Kalman filter estimate 

Pp(o) = P(t) (80b) 

where P(t) is the error covariance of the Kalman filter at the 
"present" time t. 

The two models developed for the vertical and horizonatl motions 
have been used to study the predictability of the ship motions. 
Figures 5.1 through 5.5 show simulation results assuming perfect 
state knowledge. Similarly perfect state knowledge has been 
assumed and the covariance has been propagated using equations (79) , 
(80). Figure 5.6 is a plot of the heave, pitch motions rms error 
versus prediction time, while Figure 5.7 depicts the sway, roll. 



93 


yaw rms erros versus prediction time. The error has been non— 
dimensionalized with respect to the corresponding rms motions. 

As expected, the error tends to 100% for large prediction 
times. Roll is a narrow band process and as expected it is the 
most predictable motion, up to ten seconds ahead. The remaining 
four motions are predictable up to five seconds ahead. 

These results hold in the ideal case. The actual performance 
will be lower due to the presence of noise in the measurements, 
fewer measurements than states and modeling errors. 

To assess the effect of measurement noise some simulations were 
made, whose results are shown in Figures 5.8 through 5.10. A rather 
extreme case was considered: In the case of the vertical motions 

only two measurements were available, heave and pitch, and in the 
case of the horizontal motions, only the sway, roll, and yaw motions 
were available. As seen in the figures, the same noise used to 
derive the Kalman filter gains was used, which is quite significant. 
As expected, the performance deteriorated although roll is still 
predictable up to eight seconds. The other four motions are pre- 
dictable up to about two seconds. 

In the case of modeling errors, let the correct model be 

X = AX + W (81) 


while the prediction model is 


X = A* x^ 
-p -P 


(82) 


with A* = A + 5 A. Then the error e = x, 
equation 

e = A*e + 5 A X - W 


X is governed by the 

(83) 
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or by letting = 


X 

k *0 


^1 



' A* 

5 A 


' -W 

rr 1 

0 

L* 

A 

+ 

W 


(84) 


so by denoting: 

P = E{ee^}, V = E{xe'^} ,V = E{xx'^} (85) 
the following equations are obtained 


O 

P 


o 

V 


o 

u 


A* P + P A*^ + 5A*V + v'’^*6A^ + V, 


A V + V A*'^ + U*5 a"^ - V, 


= A*U + U A + V, 


(86abc) 


V U, the covariance of the vessel model, can be assumed to be in 
steady state so the first two equations can be used to propagate 
V and P. Figures 5.12, 5.13, and 5.14 depict the error covariance 
of the vessel motions when the model used is different than the 
actual one. Since one of the most critical parameters is the modal 
frequency, its effect has been studied: The nominal value is 0.52 

rad/sec, while the value used in the prediction filter is 0.72 rad/ 
sec. The covariance at the initial time is assumed to be zero 
(perfect state knowledge) . 

As can be seen from these figures, the effect of modeling error 
is important in the case of the modal frequency, providing a re- 
duction of about 30-50% in the prediction time within prescribed 


confidence limits. 
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The sinusoidal behavior of the covariance propagation at 
about twice the motion natural frequency (as seen in the case of 
roll, for example) can be explained by the form of the covariance 
equation 

P = AP + Pa'^-V (87) 


For example the unforced equation 


O rp 

P = A P + P A 


P(0) = Q 

has the solution 

T, At„ A^^t 

P (t) = e Qe 

which is composed of exponentials in the form 
^(^i + ^j)t 


(88ab) 


(89) 


where are the eigenvalues of A. As a result 2X^|^ will appear and 
in the case of roll it is obvious that twice the roll natural fre- 
quency dominates the response. 
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Reduced Number of States 

Another aspect of interest is omitting states. In such a case 
we denote by x the nominal state of dimension n and by x* the 
implemented state of dimension m, r^m. We assume that x* is 
obtained by simply omitting some states, so that 

X* = Wx (90) 


where 


w 



(91) 


with I the unit (mxm) matrix, 
m 

X = A X + BW^ 

£ = Cx + W 2 


The nominal system equation is 

(92ab) 


while the prediction filter is 


x^= A* X 

-p -p 


% 


= ^p 


(93ab) 


where x has dimension m. A* is the mxm reduced system matrix and 
P 

C* the reduced observation matrix. Then we define; 

5a = w"^ A* W -A (94) 


and we obtain the covariance equations [8] : 

P = w’^«A**W«P + P»W^»A*'^«W + 5A«V + v'^KSa'^ + 

V = A*V + V«w'^»A*'^«W + U^Sa"^ - 
° = AU + UA*^ + 


(95abc) 
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where P, V, and U are the same quantities defined before, eX' 
cept that the error is defined as; 


T 

e = W X - X 
-P 


(96) 


Note that we may eliminate any row from the original system, 
so its is convenient to set all the rows to be eliminated at the 
end by performing a row permutation. The matrix that interchanges 
rows has the form 


S 


00 ... 010...0 


(97) 


where each row has only one nonzero entry equal to one, and no row 
is the same as any other one. 

The inverse of S is thus minimizing the computational effort 
Once the appropriate permutations have been established, it is easy 
to construct S and equation 


X = Ax + BW^ 


(98ab) 


becomes 


where 


z = A Z + B W, 
- P - P-1 

y = C z 
^ P - 


z = Sx 

T 

A = S A S 
P 

B = S B 
P 



(99ab) 


(lOOabcd) 



98 


Then we proceed to determine the error covariance as explained 

above using A , B , and C as the system matrices. 

P P P 

The inclusion of non-minimum phase zeros was considered to be 
an important part of the overall modeling. This was confirmed by 
studying the effect of omitting these zeros on the prediction error 
covariance. This is seen in Figure 5.15, where the heave and 
pitch rms error, non-dimensionalized over the corresponsing rms 
motion, is plotted versus prediction time. As expected, pitch error 
increases substantially, since pitch lags heave at low frequencies 
by 90°. Because heave and pitch are coupled, the error in heave 
is also affected, resulting in poor prediction of both heave and 
pitch. 
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Figure 5.1 

Heave Simulation Results and Its Prediction (dotted line starting 
at t=40 sec) for U=21 ft/sec and <J)=0® , and in sea state 5. 

Perfect state knowledge is assximed. 
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Figure 5.2 

Pitch Simulation Results and Its Prediction. Same conditions 

as in 5.1. 
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Figure 5.3 

Sway Simulation Results and Its Prediction (dotted line starting 
at t=40 sec) for U=15.5 ft/sec and (j)=45® , and in sea state 5. 
Perfect state knowledge is assumed. 
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Figure 5.4 

Roll Simulation Results and Its Prediction. Same conditions 


as in 5.3. 
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Figure 5 . 5 

Yaw Simulation Results and Its Prediction. 

as in 5.3. 


Same conditions 
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Fig ure 5.6 

RMS Prediction Error Over RMS Motion Versus Pred: 
for Heave and Pitch, U=21 ft/sec, W =0.72 rad/sec 
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HEAVE SIMULATION. PREDICTION 
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Figure 5 . 8 

Heave Simulation Results, Its Kalman Filter Estimate (up to 
40 sec) and Its Prediction Using the Kalman Filter Estimate 
(after t=40 sec). Same conditions as in 5.1. 
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Figure 5.9 

Sway Simulation Results, Its Kalman Filter Estimate (up to 40 
sec) and Its Prediction Using the Kalman Filter Estimate 
(after t=40 sec). Sam conditions as in 5.3. 
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Figure 5.11 ; 

Yaw Simulation and Prediction. Same conditions as 5.9 
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Figure 5.12 

rror Versus RMS Motion Versus Prediction Time 

1 W =0.52 rad/sec, used V7 =0.72 rad/sec. All 
m m 

as in 5.3. 
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Figure 5.13 

RMS Error Versus Prediction Time Roll. Same condition as 

in 5.12. 
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CONCLUSIONS 


A satisfactory approximation of the ship motion equations 
as provided by hydrodynamic theory has been achieved. The 
approximation is valid within the wave frequency range and for 
seas described by the Bretschneider spectrum, whose major 
limitations are 

(a) uni-directional seas 

(b) unlimited fetch, deep water 

The resulting two groups of motions, i.e., heave-pitch and 
sway-roll-yaw can be approximated separately requiring 15 and 16 
states respectively. If both must be used a 25 state system is 

required. 

The model depends parametrically on the ship speed, the 
wave angle and the significant wave height and modal frequency. 

The Kalman filter is designed using as measurement noise 
intensity, values from ship vibration amplitudes. For sway- 
roll-yaw the vibration levels are small, nonetheless, to bound 
the filter eigenvalues below 2.0 rad/sec similar values were 
used, as for heave-pitch. 

It should be remembered that heave and pitch are related 
by a non-minimum phase transfer function, resulting in reduced 
filter accuracy. Actually, heave is 90° out of phase for low 
frequencies with respect to all the other motions. 
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A sensitivity analysis of the filter performance indicates 
that the most critical parameter is the spectrum modal fre- 
quency. It should be remembered that a sea spectrum may contain 
more than one peak, in which case it is essential to obtain an 
accurate estimate of both peak f requenices . 

Of particular interest is the fact that the wave direction 
does not have a significant influence on the estimation error. 

This means that although our modeling used a uni-directional 
Bretschneider spectrum, it can be applied in its present form 
for directional seas. 

The models derived herein can be used to predict the ship's 
motions up to 5 seconds ahead in time for all motions and 10 
seconds for roll. When modeling errors and noise are taken into 
account, a more realistic estimate of 2-3 seconds for all motions 
and 6-8 seconds for roll . is obtained. 

Again, the modal frequency of the sea spectrum is the most 
critical parameter. Also, the nonminimum phase zeros can deterio- 
rate the performance of the predictor significantly if omitted. 
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APPENDIX 1 
Hydrodynamic Theory 


In the text, the simple one dimensional equation of heave 
motion was derived to demonstrate the principles involved. Here 
we will proceed to write the overall equations of motion. 

We will avoid extensive hydrodynamic theory developments 
since [l^ ] and [ 17 ] provide an in depth coverage. Within linear 
theory, we intend to write the added mass, damping and exciting 
force terms. 

The equations of motion can be written as 


{-0)^ [M + a] + iu) B + C} X = Fn 


( 1 ) 


where 


M = {M,j) , A . {A..}, B = {B..}. C = {C,^} 


We can find the various matrices from hydrodynamic theory 
[15 ], [ 17] and dynamics. We omit surge as a second order 
quantity so that x is a vector of dimension 


M = 




m 

0 

-mz 


0 


O 

m 

0 


c 

0 0 


0 


-mz 0 
c 


0 

1 


0 

0 


XX 

0 I 


yy 


I 0 
zx 


0 
0 

1 


xz 


zz 


(3) 


with m the mass of the ship, z the distance of the center of 

c 
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gravity, vertically, from the origin. 


A = 


2 

o 

A4 2 
0 

As 2 


0 

As 3 
0 

As3 

0 


A 2 ^ 
0 

A 4 4 
0 

Ae 4 


O A2 6 
As 5 0 

0 A 4 6 
As 5 0 


0 


A 


6 6 


(4) 


B = 


02 2 
I 0 

042 
0 

0ez 


O 

03 3 
0 

0 5 3 
0 


02 4 
0 

04 4 
0 

06 4 


0 

B 3 5 
0 

0 5 5 

O 


02 6 
0 

04 6 
0 

06 6 


(5) 


C = 


0 pgA 


up 


0 

0 


0 

C 3 5 


0 A(GM)^ 0 


0 

0 

0 


0 Css 0 A(GM) 0 

L 


0 


( 6 ) 


With A the waterplane area of the ship, A the displacement, 
up 

(GM) the transverse metacentric height and (GM) the longitudinal 

T L 

metacentric height. 

The following relations hold true within strip theory [17] 
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A 33 

= 

A° 

3 3 

— 

U , A 

B 3 3 

= 

3 3 

+ 

rT A 

U a 3 3 

A 3 5 

= 

A“ 

3 5 

- 

U ^ U ^ „A 

""a 

B 3 5 

= 

B" 

3 s 

+ 

0 a - A 

UA 3 3 - Ux^ a 3 3 - b 3 3 

A 53 

= 

A® 

3 5 

+ 

“ b”,, + ^ X, b*3 
0) ^ ^ w A 

B S 3 

= 

B“ 

3 5 

- 

UA 33 - Ux^ a^3 

A 5 5 

= 

0 

A 5 5 

+ 


B 5 5 

= 

0 

B 5 5 

+ 

“ 2 A . , A 

^ B 33 + ux^ «3 3 + ;;p^A 


The superscript 0 denotes quantities at zero speed. is 

A 

the distance to the aftermost cross-section of the ship and j ' 

. are its sectional added mass characteristics. These last two 

quantities are important only for cruiser stern ships . The A? ^ , 

B?j can be found using the M.I.T. Seakeeping Program [27]* 

Similarly, we can find: 

TV Tv° U .A 

A22 " A22 “ 2 

o ^ 

B22 ~ B2 2 + U CX2 2 
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A2 4 

= 

A4Z 

= 

Az 4 - 77 Dz 4 
U) 

Bz 4 

= 

B42 

= 

0 

Bz 4 + U ciz 4 

Az 6 

= 

0 

Az 6 

+ 

U ® U ..A , A 

-2 B22 - b22 + 

Bz 6 


0 

Bz 6 

- 

0 A A 

UAz 2 + Ux^ Ctz 2 + ^ t»2 2 

1 

A4 4 

= 

0 

A4 4 

- 

U , A 
E? °2 4 

B4 4 

= 

0 

B4 4 

+ 

A * 

U tt4 4 + B4 4 

•fg 

with Bit 4 the equivalent 

nonlinear damping 

A4 6 


0 

A4 s 

+ 

u U 1.A , ^A 

■p. Bz 4 - p bz 4 + p az 4 

Bit 6 

= 

0 

B4 6 

- 

0 A A 

UAz 4 + ^^A ^ 0 ? 

Ae 2 

r= 

0 

Az 6 

- 

^ Bzz - p bzz 

Bsz 

= 

0 

Bz 6 

+ 

UAz 2 + Ux^ af 2 

As 4 

s 

0 

A4 6 

- 

U U uA 

-2 Bz 4 - -2 x^ bz 4 

B$ 4 


0 

B4 6 

+ 

UA2 1, + Ux^ af'4 
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‘■6 6 


O ^ U 2 kA ^ A 

A & 6 ~2 ^22 ~ ~2 X , 02 j + ~2 ^12 

0 5 A td A 


( 1 ) 


B 


6 6 


o i]2 o I A A 

Bee B22 + Ux^ 0.22 + ^22 


F2 = «o P / (fa + hz) <35 + a„ p hf 


T ? 

^ 6 


«o P / (f 2 + ha) + ha ] dc + «o P 3 ^ hf 


F 4 = Oo P / (f<* + h^) dC + tto P 


U 


10 ) 


F3 = P tto / (£3 + ha) dC + P a„ ^ 


^ ^A 
10 ) 


Fs = -P a, / [?(f 3 + h 3 ) + h 3 ] - P a„ 3^ x^ht 


where f ^ is the sectional Froude-Kryloff force and h^ the sectional 
diffraction force. 



4.6 ft. 


Vertical Center of Gravity = - 

(from waterplane) 

Metacentric Height = 4.16 ft. 

Longitudinal Center of Gravity = 1.07 ft. AFT 
( from admidships) 

Displacement = 6,800 ton 


Block Coefficient 


0.461 
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APPENDIX 2 

The M.I.T. Sea-keeping Program [27] was used to 
derive the hydrodynamic data. The values used to develop 
the simple models for heave-pitch and sway-roll-yaw (as 
already mentioned, to first order heave and pitch are 
uncoupled from sway, roll and yaw). 

All units are consistent such that the forces 
are obtained in tons, the moments in ton* ft., the 
linear motions in ft. and the angular motions in radians. 
The programs change the angular motions and express them 
in degrees only in the final (output) stage. 

Heave-Pitch Characte ristics 

587 
260 
15500 

4.20*10^ 

9.53*10^ 

A1 = 550 

A2 =120000 

The principal ship characteristics have as follows 

Lgp = Length between perpendiculars = 529 ft. 

B = Beam = 55 ft. 

T = Draft = 18 ft. 


M 

214 

*^33 

A°33 = 

281 

®33 

A°35= 

15500 

®35 

S5 = 

260 


T = 

“55 

3.76*10^ 

^55 

®55 

3.8*10^ 

S5 
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Sway - Roll - Yaw Characteristics 


^4 


22,800 



^44 

= 

104,000 



®44 

= 

800 + B * 
44 



C44 

= 

28,800 



^24 

= 

-760 

^24 = 

-50 


= 

181,000 

®46 

5,600 

^22 

= 

220 

®22 ^ 

10 

^6 

= 

4 . 16 * 10 ^ 



^66 

= 

3.8 * 10 ^ 



®66 

= 

130,000 



^26 

= 

14,500 

^26 = 

370 

^2 

= 

380 



^4 

= 

2,400 



^6 

= 

23,000 
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APPENDIX 3 
Simulation 


The simulation of a continuous system 
I 

X = Ax + BW^ 


y = Cx + W 2 

where the white noise signals have intensity , 

V 2 respectively, is performed by constructing the equivalent 
discrete system 

x(T+6t) = Aj^ x(t) + 


y(t) = C x(t) + W. 


k2 


where A^^ = = I + A* 6t + A^-6t^/2 I + . . . 

„ ,t+6t A(t-T)„„ , V , 

B, W, , = / e BW, ( T ) d 

k kl ^ 1 


An approximation would be 

x(t+6t) = (I+A* 6t)x(t) + (B* 6t)Wj^^ 


y(t) = Cx(t) + Wj ^2 


where W^^2 

respectively. 


discrete white noise of intensity 
given by: 


V, 


k2 
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\2 = 

If a random number generator is provided with a range be- 
tween 0 and 1, the following relation provides (and simi- 
larly 


W. 


kl 


(RND-0.5) • 


where RND is a random number between 0 and 1. 

Higher order approximations to Aj^ may be necessary in some 
cases. Problems appear in particular with lightly damped systems, 
i.e., when the matrix A possesses eigenvalues -a + ib where 
0<a<<b. Such for example, is the case of roll whose damping 
is typically around 10% of the critical value. 

The approximation: 


^ (-a+ib) 6t 


(l-a*6t) + ib*6t 


is valid provided 
a«& t<<l 
b»5. 

The next order term is: 

- ^.,,2 6t^ 

(-a+ib) -j- 

so for accuracy we require, 

I , 2,2, 6t^ I 

1 (a -b ) —2 I <<a6t 
ab6t^<<b6t 
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which results in a single additional condition 

2a I 


6t<< 


a -b 


It is easy to see that when a<<b then 


6t« min [ji , 2 . ] 


This requirement may be very - demanding for simulation by 
imposing an extremely small time step. If an appropriate step 
is not chosen then the real part of the (neglected) second order 
term reduces the first order real part, thus resulting in a re- 
duction of the already small damping. By using a higher order 
approximation, in the form 


A. (St 
e 


I + A- 6t + 


Aiugtm 


ml 


such problems are resolved. In the present case, the sway- 
roll-yaw model requires such treatment for efficient simulation. 
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Appendix 4 


Computer Program Listing 



FILE: ^■OL.EL 


FORTRAN A 


V'VI/SP CC^'.ER'.ATIOl^lAl. MONITOR SY?IF-M 


PAGE 001 


DIMENSION A( 15,15) ,C(2,1 5) 
FX(X)=5.*(-X»*4+1 .)*(1.+2 
WRl TE( 6,2000 ) 

I FORMAT (3X, ' I NP'JT THE SPEED 
REAl)(5 ,*)U,T1 
T=T1*3 . 14159/180. 

3=^32. 2 

WRITE (6, 2010) 

1 FORMAT (3X, ' INPUT THE SIGN. 
lEC) ' ) 

READ( 5 , * )7..Z , OMM 
ALF A= J/GiO,.liVl»COS(T) 

L1 1 = . 2 
U2 = 5. 

AA=FX( 01 ) 

J 0M=-( J1 +02) » . 5 

IF( ABS (01-U2 ) . LT. . 00C1 ♦07.) 
CC^F <( UM) 

IF( AA*CC.GT.O. ) GO TC 902 
U2= UM 
GO TO 903 

2 U'=UM 
AA-CC 


, SIGMA( 15 , 15) ,XI (2,2) MOD00010 

.*X* ALFAl-2.»ALFA*X'‘*5 M0000020 

M0D00030 

(FT/SEC) AND WAVE ANGLE (DEGREES)') MOD00040 

MOD00050 

MOD00060 

MCD00070 

MODOOOQO 

WAVE HEIGHT (FT) AND MODAL FREQ. ( RAD/SMOD00090 

MOD00100 
MOD001 10 
MOD001 20 
MOD001 30 
MOD00140 
M0D001 50 
MOD00160 

GO TO 901 MOD00170 

M0D00180 
M0D00190 
MOD00200 
MOD0021 0 
MOD00220 
M0D00230 


GO TO 903 

601 ZS-UM+UM* *2* ALFA 
GAM=ZS*2 . »♦ . 25 

ZS1 =( 1 . + ( ZS/GAM)**4) **3/2S**4 

Z5."> = EX P(-1 . 2 5/UM» + 4)/UM* *5/( 1 . +2. +ALFA + UM) 

BET =ZS 1 *2S2 

A33=495 . 

830=260. 

C33=537. 

A35=1550C. 

835=1 1 700. +280. "U 
C35=260. *U+1 7250. 

AF')3.= 15 500. 

853=1 I 700. -280. *U 
Cr>3 = -260< U + 1 7250. 

A55=7. 96< 10. **S 
B5:> = 3. 8*10*“S 

C55 = -2 30 . *U* *2 + 9.53*10.* *6 

D=A33+ A55-A35 + A53 

Dl =A55/D 

D2=-A35/D 

D3=-A53/D 

D4=A33/D 

Z1=-(D1*P33+D2*B53) 

Z2=-(D1*C33+D2*C53) 

Z3 = -(D1'‘B35 + D2*B55) 

- ( 0 1 *C35+D2 + C55 ) 

Z5 = - ( 0 3* B3.3-I 04* 853 ) 

Z6=-(D3*C33+D4*C53) 

Z7=.-( 0 3*835+04*855) 

Z8=-( D3*C35+D4*C55 ) 

PL=529. 


MOD00240 
MOD00250 
M0D00260 
MOD00270 
MOD00280 
MGD00290 
MOD00300 
MOD0031 0 
M0D00320 
MOD00330 
M0D00340 
MOD00350 
MOD00360 
M0D00370 
MOD00380 
MOD00390 
MOD00400 
M000041 0 
M0D00420 
MOD00430 
MOD00440 
M0000450 
M0D00460 
MCD00470 
M0D00480 
M0D00490 
M0D00500 
MOD0051 0 
MOD00520 
MOD00530 
MODC0540 
MOD00550 


Program to prepare 
matrices for the heave- 
pitch motions. 
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file: 


1 


800 




') 




MODEL FORTRAN A 


VM/SP CONVERSATIONAL MONITOR SYSTEM 


PAGE 002 


FG-50RT(2.*3.14159*G/(RL*CCS(T)+B))+2.»3.14159/( 

(T) 

ZT=.7071 

OMNE=QMM+OMM**2*U/32 -2*C0S{ T ) 

TH1 = (2 . +ZT«F0 + + 5+F0**4*0!VlHE) 
yH2=(rO*+4*(4,*ZT+*2“*I • )+2.*ZT +F0**3*0MNE ) 

TH3=-( OMNE*FO* *4+2. »2T*F 0**5) 

TH4=F0**4 
A '■ -550 . 

A2=120000. 

TH5=01 *41 *F0**2 
TH6=D2*A2*COS(T) 

TH7=D3*A1 *F0**2 
THS=D4*A2‘C0S(T) 

0MN=0.v1M*GAM 

S0=1 . 25/4. *ZZ**2/0MM**5*BET 
VP1 =0V,N* + 2 
VP2=0MN*2.*ZT 
VP3=VP1 +SQRT (SO) 

DO 800 1=1,15 
00 800 0=1 , 1 5 
AU.J)=0. 

A(1 ,2)=1 . 

A(3,4)=1 . 

A ( 5 , 6 ) = 1 . 

A(7,8) =1 . 

A(a,10)=1. 

A 1 1 1 , 1 2 ) = 1 . 

AM3, 1 4) = 1 . 

Af 15, 1 5)=-QMNE 
Ai8,7) =Z2 
A I 8 , 8 ) = Z 1 
Ai8,9)=Z4 
A(S,10)=Z3 
A.' '0. 7 )=Z6 
4(10,3 )=Z5 
A(i 0.9)=Z8 
A(10,10)=Z7 
A(15,1 1 )=TH1 
A( 15,1 2)=TH2 
A(15,13)=TH3 
Ar‘5,14)=TH4 
A(8, 1 1 ) =TH5 
A^8 , 1 5 ) =THS 
A(10, 11 )=TH7 
A(10,15)=TH8 
A(2.1 ) =-VP1 
A(2,2) =-VP2 
Al 2,4) =VP1 
A(4.3)=-VP1 
A(4,4) =-VP2 
A(4,6)=VP1 
A(6,5) =-VPl 
A(6,S) =-VP2 


MOD00560 

RL*CCS(T)+3)*U*COSMOD00570 
MOD00580 
M0000590 
MOD00600 
MOD00610 
MOD00620 
MOD00630 
MOD00S40 
MOD00650 
MOD00660 
MOD00670 
MCDOO 08 O 
M0000690 
MOD00700 
MOD00710 
MOD00720 
MOD00730 
MOD00740 
MCD00750 
MOD00760 
MOD00770 
MOD00780 
MOD00790 
MOD00800 
M0D00310 
MOD00820 
MOD00830 
MOD00840 
M0D00850 
M0D00860 
MOD00370 
MOD00380 
M0D00890 
MOD00900 
MOD0091 0 
MOD00920 
MOD00930 
M0D00940 
MOD00950 
MQD00960 
M0D00970 
MOD00980 
MOD00990 
M0D01 000 
MOD01010 
MOD01 020 
M0D01 030 
MOD01 040 
M0D01 050 
M0D01 060 
MOD01070 
MOD01080 
MOD01090 
MODOIIOO 
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A(12, 1 1 )=--F0**2 
A( 12, 1 2) = -2*ZT + F0 
A( 12, 1 3)^F'J**2 
A( 14. 1 3)=-F0**2 
A( 1 4 , 1 4) =-2» FO + ZT 
Al 14, 1 )=VP3 
WRITE! 6,2200 ) 

2200 FORMAT!'** A MATRIX') 

DO 890 1=1,15 

WRITE! 9,2400 ) ! A(I , J ) , J= 1 , 1 5) 
2400 F0f:MAT(eE10.4) 

E90 CONTINUE 
END 


MD001 1 1 0 
MOD01 120 
M0D01 1 30 
M0D01 140 
M0D01 150 
M0D01 160 
M0D01 170 
V0D01 180 
M0D01 190 
M0D01200 
M0D01 21 0 
M0D01220 
M0D01230 
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FILE: SRY FORTRAN A 


VM/SP CONVERSATIONAL MONITOR SYSTEM 


MODEL FOR THE SWAY ROLL YAW MOTIONS OF A 00-963 DESTROYER 
REVISED VERSION 

LAST CORRECTED SEPTEMBER 11 1981 


DOUBLE PRECISION A( 16 , 16) , XI ( 16 , 16 ) , C(3 , 16) , CV( 3 , 16 ) , CA(3 , 16) 
DOUBLE PRECISION THETA(3,3) 

DOUBLE PRECISION AH( 3 , 3 ) , BH( 3 , 3 ) 

DOUBLE PRECISION P( 3 , 3 ) , R( 3 , 3 ) , T(4 , 4 ) , U(4 , 6 ) 

DOUBLE PRECISION ZET( 16 , 16 ) , WR( 16) , WI ( 16 ) , APS( 16 . 16) 


SPEED, ANGLE AND SEA STATE INPUT 


WRITE (6,100) ^ 

100 F0RMAT(2X, 'INPUT THE SPEED (FT/SEC) AND ANGLE (DEGREES) ) 
READ(5,*)V,T1 

T1=T 1/57. 29578 

WRITE(6,101) , , 

101 F0RMAT(2X, 'INPUT THE SIGNIFICANT WAVE HEIGHT (FT) ', 

1'AND MODAL FREQUENCY (RAD/SEC)') 

READ(5, »)ZZ,OMM 


SEA MODEL 

ALFA=V/32. 2*0MM*C0S(T1 ) 

DELTA=SQRT( 1 . -2 . *ZETTA**2+SQRT(9 . -4 . *ZETTA**2+4 . *ZETTA**4 ) )/2 . 
OMN=ABS(OMM*( 1 .+ALFA)/DELTA) 

BET=( ( 1 . -DELTA»*2)*»2+4.*ZETTA**2*DELTA**2)»*3/DELTA»*4 

BET=ABS(BET*EXP(-1 .25)/( 1 .+2. *ALFA) ) 

S0=O.3125*ZZ**2/0MM*BET 

VP1=0MN**2 

VP2=0MN*ZETTA*2 

HYDRODYNAMIC DATA 
FORCES 

A2=310.*SIN(T1) 

A4=2120.*SIN(T1) 

A6=11300.»SIN(T1) 

ZET2=0.72*SIN(T1) 

ZET4=0.7*SIN(T1 ) 

ZET6=0.35*SIN(T1 ) 

0M2=(O.6+V*C0S(T1)/89.444)*SIN(T1) 

OM4=(O.76+V*C0S(T 1 )/55 . 748 ) *SIN(T 1 ) 

0M6=(O.96+V*C0S(T1 )/34 . 939 ) *SIN(T 1 ) 

ADDED MASS - DAMPING 

0MES=(0.425+V*C0S(T1)/178.27)**2 
AH( 1 , 1 )=223. 

AH(1,2)=-759. 

AH(2, 1)=-759. 

AH(2,2)=22900. 

BH(1, 1)=10.6 


SRY00010 
SRY00020 
SRY00030 
SRY00040 
SRY00050 
SRY00060 
SRY00070 
SRY00080 
SRY00090 
SRY00100 
SRY001 10 
SRY00120 
SRY00130 
SRY00140 
SRY0O150 
SRY00160 
SRY0O170 
SRY00180 
SRY00190 
SRY00200 
SRY00210 
SRY00220 
SRY00230 
SRY00240 
SRY00250 
SRY00260 
SRY00270 
SRY00280 
SRY00290 
SRY00300 
SRY00310 
SRY00320 
SRY00330 
SRY00340 
SRY00350 
SRY00360 
SRY00370 
SRY00380 
SRY00390 
SRY00400 
SRY00410 
SRY00420 
SRY00430 
SRY00440 
SRY00450 
SRY00460 
SRY00470 
SRY00480 
SRY00490 
SRY00500 
SRY00510 
SRY00520 
SRY00530 
SRY00540 
SRY00550 


Program 
for the 


PAGE 001 

to prepare matrices 
sway-roll-yaw motions 
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BH( 1 ,2)=-55.4 
BH(2, D--55.4 
BH(2,2)-887. 

AH( 1 ,3)=14600.+V*BH( 1. 1 )/OMES 
AH(3, 1 ) = 14600. -V»BH( 1, 1 )/OMES 
AH(2,3)=182000.+V*BH(1 ,2)/OMES 
AH( 3,2)= 182000. -V*BH(2, 1 )/OMES 
AH(3,3)=4. 18E6+V**2*AH( 1 , 1 )/OMES 
BH( 1 ,3)=423. -V*AH( 1,1) 

BH(3, 1)=423.+V*AH( 1, 1) 

BH(2,3)=6270.-V*AH(1 ,2) 

BH(3,2)»6270.+V»AH( 1 ,2) 

BH(3, 3 )= 144000. +V**2*BH( 1 , 1 )/OMES 
C NONLINEAR ROLL FACTOR 

FBV=3. 

BH(2,2)=BH(2,2)*FBV 

C 

C MASS AND SPRING CONSTANT 
C 

AH( 1 , 1 )=AH( 1 , 1 )+215. 

AH( 1 ,2)=AH( 1 ,2)+988. 

AH( 1 ,3)=AH( 1 ,3)-230. 

AH(2, 1 )*AH{2, 1 )+988. 

AH( 2, 2)=AH( 2,2)+ 104000. 

AH(3, 1 )=AH(3, 1 )-230. 

AH(3,3)=AH(3,3)+3.76E6 

044=28800. 

C 

C CONSTRUCTION OF THE MATRICES T AND U 
C 

CALL RINV(AH.R) 

CALL RMUL(R,BH,P) 

DO 400 1=1,4 
DO 401 J=1 ,4 
401 T(I,J)=0. 

DO 400 J=1 ,6 
400 U(I ,0)=0. 

T( 1 , 1 )=R( 1 ,2)*P(2, 1 )/R(2,2)-P( 1 , 1 ) 

T( 1 ,2) = R( 1 ,2)/R(2,2) 

T( 1 ,3)=R( 1 ,2)*P(2,2)/R(2,2)-P( 1 ,2) 

T( 1 ,4)=R( 1 ,2)*P(2,3)/R(2,2)-P( 1 ,3) 

U( 1 , 1 )=R( 1 , 1 )-R( 1 ,2)*R(2, 1 )/R(2,2) 

U( 1 ,5)=R( 1 ,3)-R( 1 ,2)*R(2,3)/R(2,2) 

T(4,1)=R(3,2)*P(2,1)/R(2,2)-P(3,1) 

T(4,2)=R(3,2)/R(2,2) 

T(4,3)=R(3,2)*P(2,2)/R(2,2)-P(3,2) 

T(4,4)=R(3,2)*P(2,3)/R(2,2)-P(3,3) 

U(4, 1)=R(3, 1)-R(3,2)*R(2, 1)/R(2,2) 
U(4,5)=R(3,3)-R(3,2)»R(2,3)/R(2,2) 

T(3,2)=1 

T(2, 1) = -P(2, 1)*T(1 . 1 )-P(2,3)*T(4, 1) 
T(2,2)=-P(2,1)*T(1,2)-P(2,3)*T(4,2)-P(2,2) 
T(2,3)=-P(2, 1 )*T( 1 ,3)-P(2,3)»T(4,3)-R(2,2)*C44 
T(2,4)=-P(2, 1 )*T( 1 ,4)-P(2,3)*T(4,4) 

U(2, 1 ) = -P(2, 1 )*U( 1 , 1 )-P(2,3)*U(4, 1 ) 


SRY00560 
SRY00570 
SRY00580 
SRY00590 
SRY0O600 
SRY00610 
SRY0O620 
SRY00630 
SRY0O640 
SRY00650 
SRY00660 
SRY0O670 
SRY00680 
SRY00690 
SRY00700 
SRY00710 
SRY00720 
SRY00730 
SRY00740 
SRY00750 
SRY00760 
SRY00770 
SRY0O780 
SRY00790 
SRY00800 
SRY00810 
SRY00820 
SRY0O83O 
SRY00840 
SRY0O850 
SRY00860 
SRY00870 
SRY0O880 
SRY00890 
SRYOO90O 
SRY00910 
SRY00920 
SRY00930 
SRY00940 
SRY00950 
SRY00960 
SRY00970 
SRY00980 
SRY00990 
SRY01000 
SRY01010 
SRY01020 
SRY01030 
SRY01040 
SRY01050 
SRY01060 
SRY01070 
SRY01080 
SRY01090 
SRY01 100 
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FILE: SRY 


U(2,2)=R(2, 1 ) 

U(2.4)=R(2.2) 

U(2,5)=-P(2, 1)*U(1 ,5)-P(2,3)*U(4,5) 
U(2,6)=R(2.3) 

CONSTRUCTION OF THE MODEL MATRICES 

DO 402 1=1,16 
DO 402 J=1 , 16 
A(I . J)=0. 

402 XI(I.J)=0. 

XI(6,6)=S0*3. 1415926 
A(1,2)=1. 

A(2, 1)=-VP1 
A(2.2)=-VP2 
A(2,4)=VP1 
A(3.4) = 1 . 

A(4.3)=-VP1 
A(4.4)=-VP2 
A(4,6)=VP1 
A(5,6) = 1 . 

A(6.5)=-VP1 
A(6,6)=-VP2 
A(7.8)=1 . 

A(8.2)=A2*0M2**2 
A(8,7)=-0M2**2 
A(8,8)=-2.*ZET2*0M2 
A(9. 10)=1 . 

A( 1O,2)=A4*0M4**2 
A( 1O,9)=-0M4**2 
A( 10. 1O)=-2.*ZET4*0M4 
A(11,12)=1. 

A( 12.2)=A6*0M6**2 

A( 12, 1 1 )=-0M6**2 

A( 12, 12)=-2.*ZET6*0M6 

DO 404 1=1,4 

IP=I+12 

DO 405 J=1 ,6 

JP=J+6 

405 A(IP.UP)=U(I , J) 

DO 406 d=1 ,4 
dP=J+12 

406 A(IP.JP)=T(I . J) 

404 CONTINUE 

C OUTPUT MATRICES 

DO 403 1=1,3 
DO 403 d=1 , 16 
C(I,d)=0. 

CV(I .d)=0. 

403 CA(I,vJ)=0. 

C(1,13)=1. 

C(2,15)=1. 

C(3, 16)=1 . 

CV(2, 14)=1 . 

DO 407 1=1,6 


SRY01 1 10 

SRY01 120 

SRY01 130 

SRY01 140 

SRY01 150 

SRY01 160 

SRY01170 

SRY01 180 

SRY01 190 

SRY01200 

SRY01210 

SRY01220 

SRY01230 

SRY01240 

SRY01250 

SRY01260 

SRY01270 

SRY01280 

SRY01290 

SRY01300 

SRY01310 

SRY01320 

SRY01330 

SRY01340 

SRY01350 

SRY01360 

SRY01370 

SRY01380 

SRY01390 

SRY01400 

SRY01410 

SRY01420 

SRY01430 

SRY01440 

SRY01450 

SRY01460 

SRY01470 

SRY01480 

SRY01490 

SRY015O0 

SRY01510 

SRY01520 

SRY01530 

SRY01540 

SRY01550 

SRY01560 

SRY01570 

SRY01580 

SRY01590 

SRY01600 

SRY01610 

SRY01620 

SRY01630 

SRY01640 

SRY01650 
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IP=I+6 

CV(1,IP)=U(1,I) 

CV(3,IP)=U(4,I) 

CA(1.IP)=-P(1.1)*U(1.I)-P(1,3)»U(4.I) 

CA(2.IP)=U(2.I) 

CA(3,IP)=-P(3,1)*U(1,I)-P(3,3)*U(4.I) 

407 CONTINUE 

CA( 1 ,8)=R( 1.1) 

CA( 1 . 10)=R( 1.2) 

CA( 1 . 12)=R( 1.3) 

CA(3.8)-R(3. 1) 

CA(3. 10)=R(3.2) 

CA(3. 12)'R(3.3) 

DO 408 1=1.4 
IP=I+12 

CV(1,IP)=T(1.I) 

CV(3.1P)=T(4.I) 

CA(1.IP)=-P(1.1)*T(1.I)-P(1.3)*T(4.I) 

CA(2.IP)=T(2.I) 

CA(3.IP) = -P(3. 1)*T(1.I)-P(3.3)*T(4.I) 

408 CONTINUE 

CA( 1 . 14)=CA( 1 . 14)-P( 1 .2) 

CA( 1 . 15)=CA( 1 . 15)-R( 1.2)*C44 
CA(3. 14)«CA(3. 14)-P(3.2) 

CA(3. 15)=CA(3. 15)-R(3.2)*C44 

OUTPUT 

GO TO 399 

IER1=0 

IER2=0 

CALL TRNATB(16.16.16.16.A.APS) 

CALL LYPCND( 16.16. 16 . APS . XI . ZET . WR . WI . IER1 . IER2) 

DO 174 1=1.16 
XI(I.I)=-XI(I.I) 

174 WRITE(9. 175)XI(I .1) 

175 F0RMAT(D1 1.5) 

399 CONTINUE 

WRITE(9.310) 

310 FORMAT( 10X. MATRIX A ) 

WRITE(9.320)( (A( I . J) .0*1 . 16) . 1=1 . 16) 

320 F0RMAT(8E14.4) 

WRITE(9.330) 

330 FORMATdOX,' MATRIX C - ') 

WRITE(9.320)((C(I.d).J=1.16).I»1.3) 

WRITEO.340) 

340 FORMATdOX.' MATRIX XI ) 

WRITE(9.320)( (XI(I . J). J«1 . 16) . 1=1 . 16) 

DO 350 1=1.3 
DO 350 d«1.3 
350 THETA( I . J)=0. 

WRITE(6.360) 

360 F0RMAT(2X. 'INPUT THE SWAY-ROLL-YAW NOISE INTENSITIES 
READ(5.*)THETA( 1 . 1 ) . THETA( 2 . 2 ) . THETA( 3 . 3 ) 
WRITE(9.370) 


SRY01660 
SRY01670 
SRY01680 
SRY01690 
SRY01700 
SRY01710 
SRY01720 
SRY01730 
SRY01740 
SRY01750 
SRY01760 
SRY01770 
SRY01780 
SRY01790 
SRY01800 
SRY01810 
SRY01820 
SRY01830 
SRY01840 
SRY01850 
SRY01860 
SRY01870 
SRY01880 
SRY01890 
SRY01900 
SRY01910 
SRY01920 
SRY01930 
SRY01940 
SRY01950 
SRY01960 
SRY01970 
SRY01980 
SRY01990 
SRY02000 
SRY02010 
SRY02020 
SRY02030 
SRY02040 
SRY02050 
SRY02060 
SRY02070 
SRY02080 
SRY02090 
SRY02100 
SRY021 10 
SRY02120 
SRY02130 
SRY02140 
SRY02150 
SRY02160 
SRY02170 
(3 NUMBERS) ' )SRY02 180 
SRY02190 
SRY02200 
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370 FORMAT( 10X, ' MATRIX THETA ') 

WRITE(9,320)((THETA(I.iJ).d=1,3).I=1,3) 

END 

SUBROUTINE RINV(A.B) 

C 3X3 MATRIX INVERSION 

DOUBLE PRECISION A( 3 , 3 ) , B( 3 . 3 ) , C( 3 , 3 ) 

C( 1. 1)=A(2,2)*A(3.3)-A(3,2)*A(2.3) 

C( 1.2)=A(2. 1)*A(3,3)-A(3. 1)*A(2,3) 
C(1.3)=A(2.1)*A(3.2)-A(3,1)*A(2,2) 

C(2, 1 )=A( 1 ,2)»A(3,3)-A(3,2)*A( 1,3) 

C(2,2)=A(1 . 1)*A(3.3)-A(3. 1 )*A( 1 ,3) 

C(2,3)=A(1, 1)*A(3.2)-A(3, 1)*A(1,2) 
C(3,1)=A(1,2)*A(2.3)-A(2.2)*A(1,3) 

C(3,2)=A(1, 1)*A(2,3)-A(2. 1)*A(1,3) 

C(3,3)=A( 1 , 1 )*A(2,2)-A(2. 1 )*A( 1,2) 

DET = A( 1 , 1 )*C( 1 , 1 )-A( 1 ,2)*C( 1 ,2)+A( 1 , 3)*C( 1 ,3) 

B(1,1)=C(1,1)/DET 

B(1,2)=-C(2,1)/DET 

B(1,3)=C(3,1)/DET 

B(2, 1)=-C(1,2)/DET 

B(2,2)=C(2,2)/DET 

B(2,3)=-C(3,2)/DET 

B(3, 1)=C(1,3)/DET 

B(3,2)=-C(2,3)/DET 

B(3,3)=C(3,3)/DET 

RETURN 

END 

SUBROUTINE RMUL(A,B,C) 

C 3X3 MATRIX MULTIPLICATION 

DOUBLE PRECISION A(3 , 3 ) , B( 3 , 3) , C(3 , 3 ) 

DO 100 1=1,3 
DO 100 d=1,3 
C(I,U)=0 
DO 100 K=1 ,3 

100 C(I, J)=C(I,d)+A(I,K)*B(K,d) 

RETURN 

END 


SRY02210 

SRY02220 

SRY02230 

SRY02240 

SRY02250 

SRY02260 

SRY02270 

SRY02280 

SRY02290 

SRY0230O 

SRY02310 

SRY02320 

SRY02330 

SRY02340 

SRY02350 

SRY02360 

SRY02370 

SRY02380 

SRY02390 

SRY02400 

SRYO2410 

SRY02420 

SRY02430 

SRY02440 

SRY02450 

SRY02460 

SRY02470 

SRY02480 

SRY02490 

SRYO250O 

SRY02510 

SRY02520 

SRY02530 

SRY02540 

SRY02550 

SRY02560 

SRY02570 
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